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SUMMARY 


A method of numerical solution of the Navier -Stokes equations for the flow about 
arbitrary airfoils or other bodies is presented. This method untilizes a numerically 
generated curvilinear coordinate system having a coordinate line coincident with the 
body contour. Streamlines, velocity profiles, and pressure and force coefficients for 
several airfoils and an arbitrary rock are given. Potential flow solutions are also pre- 
sented. The procedure is also capable of treating multiple -element airfoils, and poten- 
tial flow results are presented therefor. 

INTRODUCTION 

It is imperative in numerical solution of the Navier-Stokes equations that the 
boundary conditions be represented accurately in the finite -difference formulation, for 
the region in the immediate vicinity of solid surfaces is generally dominant in determin- 
ing the character of the flow. The pressure and forces on solid bodies are directly 
dependent on the large gradients that prevail in this region near the surface, and accurate 
pressure and force coefficients require that these large gradients be represented accu- 
rately. This problem is accentuated at higher Reynolds numbers as the gradients become 
more severe. 

Therefore, almost all numerical solutions of the Navier -Stokes equations generated 
to date have treated bodies for which a natural coordinate system is available - circles, 
ellipses, spheres, Joukowski airfoils, and so forth. (Natural coordinate systems as 
defined here are those for which the body contovir under consideration coincides with a 
constant coordinate line.) The paper by Mehta and Lavan (ref. 1) has given a solution 
about a modified Joukowski airfoil accomplished by generating a natural coordinate 
system with a conformal Joukowski transformation and solving the Navier-Stokes equa- 
tions on this system. The basic Joukowski transformation was modified somewhat by 
roundii^ the trailing edge and contracting the coordinates near the body. Only one case 
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was run — a stalled flow at a 15° angle of attack and a Reynolds number of 1000. The 
method is limited to those bodies which can be generated by the Joukowski transformation 
(symmetric and cambered Joukowski airfoils, flat plates, and circular and elliptic cylin- 
ders) and does not have general applicability. Arbitrary two-dimensional bodies have 
not been successfxilly attacked as yet, primarily because of the difficulty of accurate rep- 
resentation of the boimdary conditions and the large gradients near solid surfaces when 
no coordinate line is coincident with the body contour. Some sblutidns have been •' • 
atteihpted with' interpolation between gridpbihts for boundaries not coincident -with 'coor- 
dinate lines, but this neces'sarily introduces irregularity into an otherwise smooth bound- 
ary arid places the most inaccurate difference representation in precisely the region of 
greatest sensitivity. Dawson (ref. 2) attempted to create a method for general bodies by 
the use of two uniform rectangular grids: -a fine, inner grid surrounding- the. body and 
extending for perhaps .one characteristic body dimension, and a coarse outer. grid sur-:.j. 
rounding .the inner grid and extending outwrird for perhaps 10 to 12 body, diameters-. The 
twOigrids overlap to allow for accurate transition between the , two mesh systems., ,Qnly . 
a clrcular cylinder solution was attempted, and this solution was restricted to small ; 
Reynolds numbers (R ^ 1000) because of boundary instabilities. .. - ' ' 

A method of automatic numerical generation of a' general curvilinear coordinate 
system with coordinate lines coincident with all boundaries of a generarmulticdnnected 
region contaihirig any number of arbitrarily shaped bodies has, however, been developed 
which should alleviate this problem with arbitrary bodies (ref. 3). The curvilinear coor- 
dinates are generated as the solution of two elliptic partial differential equations with 
Dlrichlet boundary conditions, one coordinate being specified to be constant oil each of 
the boundaries, and a distribution of the other being specified along' the boundaries. 

These equations are solved in finite -difference approximation by successive over- 
relaxation (SOR) Iteration. No restrictions are placed on the shape of the boundaries, 
which may even be time dependent, and the method is not restricted to two dimensions or 
single bodies. Coordinate lines may be concentrated as desired aloi^ the boundaries, 
facing of the coordinate lines encircling the body may be controlled by adjusting param- 
eters in the partial differential equations for the coordinates. 

Regardless of the shape and number of the bodies and regardless of the spacing of 
the curvilinear coordinate lines, all numerical computations, both to generate the coor- 
dinate system and subsequently to solve the Navier -Stokes equations on the coordinate 
system, are done on a rectangular grid with a square mesh, that is, in the transformed 
plane. It is also possible to cause the natural coordinate system to change in time as 
desired and still have all computation done on the fixed rectangular grid with square 
mesh. This allows the curvilinear coordinate system in the physical plane to deforni 
with a deforming body, blast front, shock, free surface, or any other boundary, keeping a 
coordinate line always coincident with the boundary at all times. The physical coordinate 
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system has been, in effect, eliminated from the problem, at the expense of adding two 
elliptic equations to the original system. 

Since the curvilinear coordinate- system has coordinate lines coincident with the 
.surface contours of all bodies present, all boundary conditions may be expressed at grid 
points. Also, normal derivatives on the bodies may be represented by using only finite 
differences.between gr^ points on coordinate lines, without need of , any interpolation, 
even though the coordinate systeni is not orthogon^ at the boundary. Numerical solu- . . 
tions for the lifting and nonlifting potential flow about Kirmlh-Trefftz airfoils obtained 
with this coordinate -system generation show excellent comparison with, the analytic. . ^ 

solutions. . . . 

This method of automatic body-fitted curvilinear coordinate generation' has been'. . - 
used to construct a finlte-^hfference solution of the fully incompressible, time -dependent; 
Navier -Stokes equations for the laminar viscous flow about arbitrary two-dimensional' '" 
airfoils or any other two-dimensional body (ref. 4). The Navier -Stokes equations are- 
written in the vorticity— stream -function formulation, with the vorticity bn the body being 
determined by a type of false -position iteration so that the no -slip boundary condition is 
satisfied. 'The solution is implicit in time, the vorticity, and the stream -function equa- 
tions being solyed simultaneously at each time step by SOR iteration. A method of con- 
trolling the spacing of the coordinate lines encircling the body has been developed in 
order to treat high Reynolds number flow, since the coordinate lines must concentrate 
near the surface to a greater degree as the Reynolds number increases. The solution is 
designed to provide the velocity field, the surface -pressure distribution, and the lift, 
drag, and moment coefficients. Results are given for separated flow over two airfoils 
and an arbitrary rock. Initial application to multiple airfoils has also been made. 

SYMBOLS 

a,b,c,d coefficients in equations (5) 
axial -force coefficient 
Cj) drag coefficient 

Cjjp friction -drag coefficient 

Cop . pressure -drag coefficient 

Cl lift coefficient 
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Cn 

Cp 

c* 

D 

D* 


ds 

E 

F 

f 

l, j 
U 

J 

k 

M,N 

m, n 
S 


P,Q 


P 

R 

S 


normal-force coefficient 
pressure coefficient 

pressure coefficient referenced to trailing -edge pressure 
differential operator; two-dimensional region (fig. 1) 
rectangular region (fig. 1) 

increment of arc tength along body surface . ^ 

maximum norm 
force on body 
function 

computational grid points; i » 1 . . . I; j » 1 . . . J 

vuiit vectors 

Jacobian 

iteration coxmter 

summation limits (eqs. (5)) 

indices 

vmit vector normal to body surface 

amplitude factors (eqs. (5)) 

pressure 

Reynolds number 

body surface 


472 


stress vector on body surface 


Tn 

t 


tn 


y 


x,y 


time 

current time 


velocity 


tangential velocity component ■ 

ph 3 TSical coordinates (nondimensionklized by'airtoil chord) 


aAY 


^1, Tg, . 


6 

0 

\ 

i,v 

a,T 

Ct> 


coefficients of natural coordinate transformation (eqs. (3)) 
Fg curves in physical plane 

curves in transformed plane 
convergence factor 
relaxation factor 
free -stream angle of attack 
coefficient in stream -function equation 
transformed coordinates 
coefficients in equation (9a) 
stream function 
value of at body 
vorticity 

value of 0 ) at body 
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Superscripts: 

L lower surface 

U upper surface 


transformation 



Subscripts: 


L V' lower surface ' 

...... • - ^ ^ K • ' 

o ' trailing edge (fig.4) 

jCJ. ..?i . V. \ 

T.Ei'' ‘ trailing edge ■ 

U • -“ upper surface - • 






1? >■’ 


x,y • . differentiation with respect to ;x -or ”y 







^,r) differentiation with respect to | or 7 ; 

oo free stream 


BODY-FITTED CURVILINEAR COORDINATE SYSTEM 


Mathematical Development 

Let it be desired to transform the two-dimensional, doubly connected region D 
boiinded by two closed contours of arbitrary shape into a rectangular region D*, as 
shown in figure 1. The general transformation from the physical plane [x,j 7 ] to the 
transformed plane [ 4 , 77 ] is given by | = ^(x,y), rj = 77(x,y). Similarly, the inverse 
transformation is given by x = y = y{i,t)). Derivatives are transformed as 

follows: 


. _ Mf,y)/3(^,T7) -y|t7? 

* a(x,y)/a(|,7/) ■ j 


(la) 


, _ 8(x,f)/a(|,7;) _ + xgf^ 

y a(x,y)/a(|,77) J 


where J is the Jacobian of the transformation J = X|y,j - Xj^y^. 
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Since the basic idea of the transformation is to generate transformation functions 
such that all boundaries are coincident with coordinate lines, the natural coordinates 
[{.il are taken as solutions of some suitable elliptic boundary value problem with one 
c/t tbese coordinates constant on the boundaries. Using Laplace’s, equation as the gener- 
ating elliptic system gives 

(2a) 

%x + ^yy ~ ® (2b) 

with Dirichlet boundary conditions: 17 = Constant = 77 j on F j , 77 = Constant = 772 on 
F2, and |(x,y) a multiple -valued solution with a branch of |(x,y) specified (but not 
constant) on F^ and F2. The curve F]^ on the physical pl^e trwsforms to the 
lower boimdary F^ of the transformed plane. Similarly, F2 transforms to F^, and 
so forth. The right and left boundaries of the rectangular transformed plane F^ and 
r| are coincident in the physical plane. The curve which transforms to these bound- 
aries connects F^ and F2 and determines a branch cut for the multiple. -valued func- 
tion |(x,y). Thus the functions and all derivatives are continuous across this cut. 

Now since it is desirable to do all numerical computation in the rectangular trans- 
formed plane, it is necessary to interchange the dependent and Independent variables in 


equations (2). Thus 

ax|| - 2/3x|^ + yx^ = 0 (3a) 

oey|| -20y|^ + yy^ = O (3b) 

where 

Of = Xt,2 + y^2 (3c) 

0 = X|X,^ + y^y^ (3d) 

y = X|2 + y^2 (3g) 


with the transformed boundary conditions: x = fiUtVi) on Fj, y = gi(|,?7i) on Fj,, 

X = o** ^2» y ” S2(^>’^2) ^2* (^ present application, x and y 

are nondimensionaliztjd with respect to the aiHoil chord.) 

The natural coordinate system so generated has a constant 77 -line coincident with 
each boundary in the physical plane. The ^ -lines may be spaced in any manner desired 
around the boundaries by specification of |x,^ at the equispaced 4 -points on the 77^- and 
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?72 -lines of the transformed plane. Control of the spacing of the t; - lines may be exer- 
cised by varying the elliptic system of which | and 77 are solutions. 

Extension to Multiple Bodies 

The same procedure for natural coordinate generation may be extended to regions 
that are more than doubly connected, that is, have more than two closed boundaries or, 
equivalently, more than one body or hole within a single outer boundary. The transfor- 
mation to the rectangular field Is illustrated in figure 2. 

The method requires that the 77 -coordinate be equal to the same constant on all the 
interior boundaries, that is, on all bodies in the field. Let all the bodies be connected 
by arbitrary cuts and, similarly, one body be connected to the outer boundary by an arbi- 
trary cut. Since the 77-coordlnate is equal to the same constant on all the bodies, it is, 

Of course, equal to that constant on the cuts between the bodies also. By contrast, the 
I -coordinate is taken constant on the cut between the body and the outer boundary. Since 
the locations of these cuts in the physical plane are not specified, the specification of 77 
or 4 as constant on a cut does not overspecify the elliptic problem^ 

Note that all bodies except one are split into two segments. Each cut appears 
twice on the transformed field boundary, the two segments, of course, corresponding to 
the two "sides" of the cut in the physical plane and thus being reentrant boundaries with 
the functions and all derivatives continuous thereon. Thus x and y have been speci- 
fied on the portions of the lower boundary of the transformed field that correspond to the 
bodies - r|^ and for the right body and F| for the left body - and also on the 
entire upper boundary, corresponding to the outer boundary in the physical field. The 
remaining portions of the lower boundary and the entire side boundaries are reentrant 
boundaries and, thus, neither require nor allow specification of [x,^ thereon. 

Again an elliptic Dirichlet problem is solved to generate the natural coordinates 
[x,^, as in the previously considered case with only a single body. All computations, 
both to generate the coordinates and subsequently to solve the partial differential system 
of interest, are again done on the rectangular field with square mesh in the transformed 
plane. 

Numerical Solution 

The relation between the transformed and physical fields for a single airfoil is 
shown in figure 3(a). The physical coordinates of I points describing the body surface 
[x,^ provide the boundary conditions along the j = 1 line; those of I points on the 
physical remote boundary, usually a circle with radius 10 or more chords, svpply the 
boundary conditions along the j = J line of the transformed field. Since the side bound- 
aries of the transformed field are reentrant, corresponding to the cut in the physical 
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plane, then = and fj+ij = ^2,} all j. Note that the values of x and y 
are not specified on these side boundaries. All derivatives in equations (3) are approx- 
imated by second-order, central -difference expressions (A^ and Arj are both unity by 
construction, the actxial values of 4 and rj being immaterial) : 

/ ^^)ij 

WiJ * 

?w)ij * ^U-1 

The resulting set of 21(J - 1) nonlinear difference equations, two for each point 
for i = 1, 2, . . .,1-1 and j = 2, 3, . . ., J - 1, were solved by accelerated Gauss- 
Seidel (SOR) iteration. The iteration was considered to have converged when the maxi- 
mum absolute change on the field between iterations was less than 10 ~5. A range of 
acceleration parameters was examined, and a value of 1.85 was nearly optimum for the 
bodies considered. 

The relation between the transformed and physical fields for two airfoils is shown 
in figure 3(b). The physical coordinates of body 2 at points i = 1 . . .11, those of 
body 1 at points i = 12 . . .13, and finally the remaining points i = 14 ... I on body 2 
are input as boundary conditions on the j = 1 line in the transformed plane. The 
remaining points i = (II + 1) . . . (12 - 1) and i = (13 + 1) . . . (14 - 1) on the j = 1 
line are reentrant points corresponding to the cut between the bodies in the physical 
plane. Therefore values at these points are not specified, but rather the relations 

^Il+k,l “ ^I4-k,l ^ll+k,0 “ ^I4-k,2 1 • • . (12 - II - 1) hold. The rest of 

the procedure is imchanged from the case of a single airfoil, except that two difference 
equations at each of the points [i,l] for i = (II -i- 1) . . . (12 - 1) are added to the sys- 
tem, so that the total number of equations is now 2I(J - 1) + 2(12 - II - 1). 

Control of Coordinate System 

Several procedures for controlling the spacing of the coordinate lines in the field 
are available and the general philosophy of such control is discussed in reference 3. 

One particularly effective procedure is to add exponential inhomogeneous terms to the 


(4a) 

♦ 

(4b) 

(4c) 

.m 

(4e) 

[‘4 
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Laplace eqiiatlons for the curvilinear coordinates, so that the coordinates are generated 
as the solutions of 



* I - in)^® + (” - ”n)2] o ' (5b) 

where the amplitudes and decay factors are not necessarily the same in the two equa- 
tions. Here the first terms have the effect of attracting | -lines to the -lines in the 
4-equati(m, and attractii^ 77 -lines to the 77,j^ -lines in the 77 -equation. The second terms 
cause I -lines to be attracted to the points [4n>’7n| ^ -equation, with a similar 

effect on 77 -lines in the 77 -equation. 

In the transformied plane these equations become 

+ YXjjrj = -J^(Px^ + (6a) 

ay^ - 2/3y|^ + yyrjr} = + 0^1]) (6b) 

POTENTIAL -FLOW SOLUTION 

Laplace Equation and Boundary Conditions 

The two-dimensional irrotational flow about any number of bodies may be described 
by the .Laplace equation for the stream function )//: . 

^xx ^yy =0 (7) 

with boundary conditions: 

On the body surface, 

Xx,y) = (8a) 
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At infinity. 


i//(x,y) = y cos 0 - x sin 0 


(8b) 


where 0 is the angle of attack of the free stream relative to the positive x-axis. Here 
the stream function is nondimensionalized relative to the airfoil chord and the free- 
stream velocity. When transformed to the curvilinear coordinate system, this equation 
becomes 




axl/^ - + Y^rjT) + ® 

where a, p, and y are defined by equations (3c) to (3e) and a and t are given by 


a = 


y|(Dx) - x^(Dy) 
1 “ 


■){ 




^ Xq(Dy) - y7;(Px) 


‘ (9b)- 
(9c) 


with 


■■ bx s ‘ . ■ ' (9d)' 

Note that Dx and Dy, and hence a and t, vanish when no coordinate contraction is 
used, that is, when the generating system' is simply equations (3). The transformed 
boimdary conditions are 


On rj = rji on 

‘ (10a) 

On 11 = ^2 ,yv V 

The uniqueness is implied by insisting that the solution be periodic in -<« < ^ < <»^ 

^ri 2 > The coefficients Qt, 0, y, cr, and f are calculated during the generation 
of the natural coordinate system. For the approximation of equations (9), second-order; 
central differences are used for all derivatives, and the resulting difference equation is 
solved by accelerated Gauss -Seidel (SOR) iteration on the rectangular transformed field. 


The solution of equations (9) on the transformed field is constructed in the same 
manner that has been previously described for the solution of equations (3). The single 
equation (9a) replaces the two equations (3a) and (3b), and the boundary conditions are 


479 



given by equations ( 10 ). The total nxunber of difference equations thus is ,I(J 7 1) for a 
single airfoil and I(J - 1) + (12 - II - 1) for twp airfoils. 


Kutta Condition 

The boundary value of ^ on the body ^ is determined by imposing the Kutta 
condition. The Kutta condition arises from physical considerations and basically asserts 
that the flow must leave the sharp trailing edge of. an airfoil section in a smooth fashion. 
In a mathematical sense this smoothness condition is guaranteed by insisting that the 
velocity on the surface of the airfoil be continuous. The continuity implies that the limit 
of tlie velocity at any point oh the 'svirface'exlsts'and is the hame regkrdlesh of the path 
along which' this point is approachetl “in partlcuiaf,' the velocity at the trailing* edge of 
the airfoil must be the same when approached from the \q)streahi ‘direction hiohg'tHe 
^ upper and lower surfaces. It is easily shown that the above ideas imply tlut the trailing 
edge is a stagnation point for airfoils having an included trailing -edge angle greater than 
zero, but only a common (possibly nonzero) upper and lower surface velocity limit is 
required for cusped trailing edges. - The common^limit condition has also been applied 
by Giesing (ref. 5) in a solution utilizing superposition of singularities. 

' Since the normal velocity component vanishes identically on the airfoil surface, 
only the tangential velocity! component need be considered.' If. ' ' is the component 

of V tai^ent to a constant 77 -llne> then ' . ! . . 

/ , . V ; y, ' 


Vt' 


( 11 ) 


On the surface the ^-derivatives are approximated by the second-order, central- 
difference expressions of equation (4a), as in the interior of the field, at all points except 
those on the cut i = 1 and i = I, where second-order, one-sided expressions are used. 
Thus 


®1,1 = i("^3,l + 4^2,1 -3f 1,1) 


( 12 a) 

( 12 b) 


The 77 -derivatives on the surface are approximated at all points by similar one-sided 
expressions: 

t 

To implement the condition of a common velocity limit numerically, the tangential veloc- 
ity component at the airfoil trailing edge is approximated by a three -point, quadratic 


480 



extrapolation in which the three points on the airfoil i^rface immediately ac^acent to the 
trailing>edge point on both xqpper and lower surfaces are utilized. This procedure is 
illustrated in figure 4. The extr^olated values are 




I j .. 4* ; - / - 

• ' ' V ► 


» - t:! !,<.• . 

(L) = 




(13a) 


(13b) 


where the subscripts, o, . lU, 2U, ,3U, IL, ,,2L,.and 3L refer to the ^-field posi- 
tlon.as indi^catj^ in f^ure,4. , .All.i/rfield position indie, es are of course unity^ The 
common-limit condition is then . . , , .. . . . . . 




(L) = 




(14) 


Siqaerposition of Solutions ' r 


Since the system to be solved is linea^r in <//, the solution for a single airfoil at 
any angle of attack may be obtained by superposing three component solutions: . (1) a 
soluUon at 0^ angle of attack with no circulation, (2) a solution at 90° angle of attack 
with no circulation, and (3) a solution with circulation but zero free -stream velocity. 
These three component solutions, written where k - 1, 2, 3, each satisfy 

equations (9), with the respective boundary conditions 

II 

o 

(i = 1 . . . I) 

(15a) 

4,j “ yi,j 

(i = 1 . . .1) 

(15b) 

o 

II 

(i = 1 ... I) 

(16a) 

II 

(1 = 1. . .1) 

(16b) 

II 

(1 = 1 ... I) 

(17a) 

o 

n 

(i = 1 . , . . I) 

(17b) 
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The complete solution with arbitrary circulation then is 

cos e + sin 6 + (18) 

The Kutta condition is then satisfied by choosing the coefficient \ such that equation (14) 
is satisfied. Thus it is only necessary to solve the system of difference equations three 
times, for a given airfoil. The solution at any angle of attack may then be obtained with- 
out re -solving the difference system. ■ 


Surface Pressure and Force Coefficients 


The pressure coefficient at' any point' in the field may be obtained from the-veloc- 
ities via4he- Bernoulli equation, which ih'the pfeserit hbhdimensional variables; is ' ' ■ ' ' 

” '• ' ■’ ■ A*' ii; I- .'i! li c. 




a ^ w 


• " 1 , ! • j • « v : M 

T ^ 


On the body surface this becomes, through use of equation (11), 


( ■; 


C - 1 2L iZ/ 2 • 

1 - j2 *^7; 


i ) ; • 


(19), 


r. ( 20 ) 


with the derivative evaluated by a second -order, one-sided difference expression. The 
nondi'mensional force on the body is* given by . 


F = -^ Cpn 


ds 


( 21 ) 


where n is the unit outward normal to the surface, and ds is an increment of arc 
length aloi^ the. surface. The lift and drag coefficients are 


CL=^Cp (-X| cos B - sin 0) d| 


(22a) 


Cd = ^ Cp(y| cos 0 + x^ sin 0)d| 


(22b) 


These integrals were evaluated with numerical quadrature by means of the trapezoidal 
rule. 


Multiple Airfoils 

With two airfoils, the boundary condition of equation (8a) is replaced by the two 
boimdary conditions: 

On the surface of body 1, 

Hx,y)^4'i (23a) 
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On the surface of body 2, 

i//(x,y) = i/zg (23b) 

Mth reference to figure 3 and the discussion in the previous section of the coordinate 


system solution/ these boundary conditions become, in the transformed field, 

i 

•^1,1 - '*'1 

K.:. ... V. M ^:i I ; I3r " ^ ^ 

‘ • i ; ' i . 

(24a) , 


(i = 1 . . .11 and i = 14 . . .1) 

(24b) 

As ins the, case of- the coordinate system solution,,the remaining portions of the J = 1 
line are reentrant bouncteries, so. that points thereon are treated as field points .rather j . 
than boundary points. The | -derivatives at the surface points 11, 12, 13, and 14 on 
the cuts between the bodies are also evaluated by using the one-sided expressions of 
equations (12) in the calculation of the velocity on the surface. 

The Kutta condition must be applied on each body. Therefore, a fourth component 
solution is added, and the four component solutions each satisfy equation (9a), with the 
boundary conditions 

o 

II 

(i = 1 . . . 11, IZ . . .13, 14 ; .,,1) 

(25a) . 

= y. T 
^i,J ^i,J 

(i = 1 . . . I) 

(25b) 

o 

II 

(i = 1 . . . il, 12 . . . 13, 14 . . . I) 

(26a) 

II 

(i = 1 . . . I) 

(26b) 

o 

II 

(i = 1 ... 11, 14 ... I) 

(27a) 

II 

(i = I2...I3) 

(27b) 


(i = 1 . . . I) 

(27c) 

II 

1-^ 

(i = 1 ... 11, 14 ... I) 

(28a) 
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(i = 12 . . . 13) 


(28b) 



= 0 (i = 1 . . . I) (28c) 

The complete solution with arbitrary circulation about each body is 

cos 9 + sin 9 + + \ 2 'P^^\i,ri) (29) 

The Kutta condition is then satisfied by choosing the coefficients A-i and \2 such that; 
equation (14) is satisfied on each-body. This requires only the simultaneous, solution of ,i , 
two linear algebraic equations. - Generalizing -to N bodies, it is necessary to solve the 
difference equation system N'+‘2 times for a given multiple airfoil system. The solu-\-, 
tion at any Orientation of the free stream may then be obtained without re -solving the 
difference system. 


Results and Comparisons 

The coordinate system for a K&rmln-Trefftz airfoil having an integral flap is 
shown in figure 5, and the streamlines and pressure distribution for this airfoil are com - . 
pared with the analytic solution (ref. 6) in figure 6. Similar excellent comparisons have 
been obtained with other K^min-Trefftz airfoils. Figure 7 shows the coordinate sys- 
tem for a Liebeck laminar airfoil, the solution for which is compared with experimen- 
tal results (ref. 7) for the pressure distribution and lift curve in figure 8. Finally, the 
coordinate system for a multiple -element airfoil is shown in figure 9, with the stream- 
lines and pressure distributions shown in figure 10. Here coordinate system control 
was employed as discussed above to attract the coordinate lines into the concave region 
formed by the intersections of the cut between the airfoils. 


APPUCATION TO THE NAVIER-STOKES EQUATIONS 


Basic Equations 

The stream -function —vorticity formulation of the two-dimensional, incompressible, 
viscous -flow equations is given by 


Wt + t//yO)x - - 


t*bcx + ^yy 
R 


(30) 


^xx + ^yy= 


( 31 ) 
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where ^ is the nondimensional stream function, the nondimensional vorticity, and 
R the Reynolds number based on the characteristic velocity and length used to nondi- 
mensionalize the basic equations. The transformed equations are 


wt + 


- 2/3w|^ + Yiarjq + 
J j2r 


(32) 


(33) 

where the coordinate system parameters a, /3, y, J, a, and r have already been 
given. Recall that these coordinate system parameters are fixed and need be calculated 
only once. 

Boundary Conditions 
- The boundary .conditions are given 
Oh the body surface, 

»//(x,y,t) = 4'\i = Constant (34a) 

|^(x,y,t) = 0 (34b) 

dn 

At infinity, 

\l/{x,y,t) = y cos 0 - x sin 0 (35a) 

o)(x,y,t) = 0 (35b) 

where n is the unit vector normal to the body surface. The function describing the 
variation of the vorticity on the body (>>b(x,y,t) is unknown and must be calculated as part 
of the solution. Initial conditions at t = 0 are those resulting from an impulsive start. 
Equations (34) and (35) may be transformed to yield boundary conditions for equations (32) 
and (33) in the transformed plane. This procedure yields the following relations: 

On = ^i.e.i on Tj), 

= ’/'b “ Constant (36a) 

= 0 ‘ (36b) 

= y(^>^2) ® "*(^>^2) (37a) 
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(37b) 


where and ri^ are the values of the 7/ -coordinate for contours rj and Fg, 
re^ectively, in the transformed plane (fig. .1). The condition specified by equation. (36b) 
gu^antees that the velocity component tangent to the transformed body surface vanishes 
on the body surface. Since the component normal to the body surface vanishes identically, 
the satisfaction of eqiration (36b) implies that the viscous no-slip condition is satisfied 
on the body surface ^i.e., along F*, or equivalently r’jJ. 

' Most solutions in the computational fluids field have relied op a modified evaluation 
of equation (33) on the boundary to determine the vorticity on the body surface 
The modification is introduced in an attempt to insure that eqviation (36b) holds — that is, 
to satisfy the no^slip condition. A variety of numerical procedures along these lines are 
documented in reference 8. The principal problem encountered with such an approach is 
that the vanishing of the tangential velocity component is implied only indirectly rather 
than. directly. . Israeli (ref, 9) has shown that these procedures are not only unreliable in 
producing a. zero tangential component, but may, in fact, even be numerically divergent. 
Israeli suggests that calculated with an iterative algorithm of the form 

(38) 

for all 1 s I S I _ i, ‘where refers to the ^-jposition along the body, tjj denotes the 
-value for contour Fj, tn is the current time, k denotes the iteration counter at 
ste^ tn, and 6 is a relaxation factor (possibly variable). Obviously, such a procedure 
can only be employed with implicit methods which require iteration of the parabolic vor- 
ticity equation at each time increment. Note that convergence of the vector sequence 

^b^\^i»7l»tn) iniplies convergence of ^^(^i>7l>^n) to that function which inherently 

satisfies the no-slip boundary condition. . ' , 

Pressure Coefficients 

If the prinxitive variable formulation of the Navier -Stokes equations (velocity- 
pressure) is evaluated on the body surface, the time derivative and Inertia terms vanish 
to yield 

yp = i v2 y (39) 

< i" ■ - V ' ' . • • , .-v- • 

where p is the nondinxensional pressure, V the.nondimensional velocity, and R ^ 
Reynolds number based on the characteristic flow parameters. Utilizing ayector iden- 
tity to eliminate y gives 
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B> = (40) 

The pressure differential in the transformed plane then is 
. dp = - yw7/)d4 . 

bitegration of equation (41) along starting at the trailing -edge position yields 

C^(€) = P(l) - Pt.E. “ 5 - yt^r/)d|* (42) 

where If is the | -value corresponding to the body trailing edge. The S]mibol Cp 
is used instead the more conventional Cp to indicate that the reference pressure is 
the trailing -edge value Pf_E. rather than the free -stream value p^^. Note that all 
quantities in equation (42) must be evaluateii bn the body sui^ace (i.e., along q » q|j. 
Central-^ference approximations. were used for all | -derivatives appearing in equa- 
tion (42), while second-order one-sided expressions were used for the 77 -derivatives. 

The numerical quadrature was performed by the trapezoidal rule. 

Force Coefficients 

. The force coefficients associated with the stress vector are obtained by integrating 
the stress vector over the body surface.. Let. F = i C;^ + j Cj^ be the total force acting 
on the body and let Tn = 1 ("^n) ^ + 1 (^0)2 stress vector on the body surface having 

outward unit normal n. Then, 

F = f TndS (43) 

'* *^S ~ 

where S is the body surface. The stress vector components (Tq^^ and (Th )2 
be expressed in terms of the primitive variables as 

(T^j = -2pn, + 4(Vi)^ H . 2j-(Vi)^ , (V2)^ja (44a) 

* ‘ ^ ’ ... ' I - 4 . 

(T5)j=-2pnj.2[(V4.(V,)^H.4(V2)ya (4ib) 

where n^[ and n2 are the x- and y-components of the normal to the body surface ii. 
The lift and drag coefficients may then be calculated by means of the conventional wind- 
axis ti^sformation as follows: 
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•The two integrals in equations (45) are referred to as the pressure and friction drag 
coefficients and are denoted C[jp and Cj^p, respectively. These integrals were eval- 
uated numerically by use of the trapezoidal rule. The 4 -derivatives were approximated 
with second -order central -difference expressions. 

Difference Equations 

A first-order backward difference is used to approximate the time derivative, 
while second-order central differences are used for the space derivatives in equa- 
tions (32) and (33)i The resulting coupled difference equations, two for each point in the 
field, were solved simultaneously by point SOR iteration at each time step. 

Implementation of the Boundary Conditions 

As indicated earlier, the basic idea used to calculate the vorticity on the body sur- 
face is to select this function so that the no -slip condition is satisfied. An 

approach suggested by Israeli (ref. 9) has already been cited in equation (38). This is 
basically the parallel -chord method (see ref. 10) and has only a linear convergence rate. 
Another method similar to false -position iteration was used to accelerate the conver- 
gence. The iterative sequence is generated by the algorithm 



where 6 is again an acceleration paranieter. The derivatives in this equation were 
approximated with second-order one-sided differences for the -derivatives and central 
differences for [the- ^-derivatives. Another method was used when numerical overflow . 
problems were encountered with the quotient in equation (46). This consisted of modi- 
f3ring the second term on the right side of equation (38) with the algebraic sign of the dlf-^ 
ference quotient. ..’;Several other approaches documented in Roache (ref. 8) ^ere also 
tried. None of these were as successful as the methods discussed above. 
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Calculation Procedure 

The vorticity and stream -function fields are converged by means of the point SOR 
technique. New boundary values of the vorticity are calculated as discussed in a previ- 
ous section. Three conditions must be met before the time step is considered to be 
converged: 


( 1 ) E{xl^,k) ^ e{xly) 

(2) E(u),k) i e(w) 


(3) 




where the terms E(»^,k), E(w,k), and E^|^,l^ are the maximum norms of the change 
E(^,k) = The terms involving e are simply the required convergence 

criteria. ^Nominal values for e(V'), e(a)), and a^re 10“®, 10"®, and 10"^, respec- 

tively.^ This procedure is repeated until convergence. Once a time step has converged, 
time is incremented and the process begins again. 


NAVIER-STOKES RESULTS 
. Solutions About Various Bodies 

To illustrate the versatility of the natural coordinate system. approach, viscous 
flows about three different bodies are presented. The bodies and associated flow condi- 
tions are 

(1) Fl^ped KirmSn-Trefftz airfoil: 0 = 15°; R = 200 

(2) Gbttlngen 625 airfoil: 0 = 5°; R = 2000 

(3) Cambered rock: 0 = 5°; R = 500 

The coordinates of these bodies are given in reference 4.. 

Several problems arose with the body vorticity calculations. At time's the iterative 
method used to calculate the body vorticity produced mildly oscillating values along the 
body boimdary. The principal cause of this result was that the method was applied point . 
by point along the boundary. Thus, the only "communication" between the body points 
was through the field iterations. This tendency was overcome in two ways: First, only 
small surface vorticity changes were allowed at each body station at each iteration. 
Second, after the new vorticity values had been calculated, a three -point weighted aver- 
age was used to smooth the new surface vorticity distribution. 
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A second problem developed with those bodies having a sharp trailii^ edge. To 
preseirve continuity of the vorticity, the vorticity at the trailing edge was held at zero^ ' ' 
Since the vorticity gradients are extremely large in the neighborhood of the trailing e^e/ 
the numerical solution had a tendency to oscillate near this point. This phenomenon is 
generally known as ''wiggles'* and, as shovm in Roache (ref. 8), is actually the solution of 
the difference equations, hi reality the wiggles are caused by the inability of the net 
function to resolve large gradients near boundaries.' / ■-. 

Flapped K&rm&n-Trefftz Airfoil 

The coordinate sy^em for the flapped K^irm&n-Trefftz airtoil profile, which pos- 
sesses a camber of 22 percent at the 0.55 -chord point, is shown in figure 5. The free- 
stream Reynolds number was taken as 200, and the flow angle of attack was 15^. Other 
data concerning the solution are given in table 1. 

Stream -function contours are given for two time steps in figure 11. The contours 
at the earlier time indicate clearly the large flow velocities over the upper suiiace of 
the airfoil and the consequent large difference in overall boundary-layer thiclmess 
between the iq>per and lower surfaces. The manner in which the zero streamline leaves 
the trailing edge indicates t^t flow separation on the upper surface is imminent. The ' 
contours for t - 1.06 illustrate a fully developed laminar separation. The boundary- 
layer thiclmess over the aft half of the upper surface has increased approximately 
300 percent. 

. In order to gain some insight into the development of laminar separation, a series 
of four, velocity prof iles are shown in figure 12. The profiles for t = 0.08 illustrate 
the iq)per -surface flow shortly after the Impulsive start. The boundary layer is very 
thin at this time. Separation has already begun at t = 0.22, as evidenced by the profiles . 
on the flap portion of the a.irfoil. Figures 12(c) and 12(d) indicate that the upper-surface 
separation, point has moved rapidly upstream to approximately the 70 -per cent -chord 
potot. Reverse flow tes been well established at t = 0.54. Note that the upper -surface 
boundary-layer thickness has increased substantially over the time span shown. 

Gbttingen 625 Airfoil 

The flow Reimolds number was 2000 at an ai^le of attack of 5^. . Additional sum- 
mary data of the solution for the Gdttingen 625 airfoil appear in table 1. The coordinate 
system sho^ in figure 13 was used in this solution. The high density of constant 77 -lines 
new the airtoil s^ace is the result of contraction to the first 15 77 -lines. In particular, 
the amplitu^ factor appearing in equation (5b) ranged from 20 000 at 77 = 1 to 13 OOO'at 
17 a 15 (increments of 500/line), while the decay factors were held constant at 1.0 for 
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17 = 1 to 14 and were 0.4 for the 15th t 7 -line. The function P(x,y) defined by equa- 
tion (5a) was set to zero as were the point -attraction parameters of equation (5b). 

. This solution developed wiggles near- the sharp trailing edge. The effect, of, the 
wiggles is dramatically illustrated in figure 14, which shows and co. conto\irs a.t . 
several times after the impulsive start. Note the distortion of the vorticity contours 
near the trailing edge. The oscillatory effects are carried upstream along the lower 
surface of the airfoil and are proceeding downward, away from the trailing edge. The 
disturbance proceeds away from the airfoil without much damping but has little effect 
on the flow in the vicinity of the airfoil after the start, as can be seen in figures 15(a) 
and (b), which show »// contours at later times. A feature of interest in figures 14(a) 
and (b) is the starting vortex which is formed and shed at the trailing edge. This vort^ 
appears just above the disturbance due to the' wiggles. 

Figure 14(c) indicates that flow separation has been initiated on the trailing-edge 
portion of the airfoil upper surface. The separation point moves rapidly upstream to 
approximately the half -chord point at t = 1.012. At this time the upstream movement 
of the separation point slows down considerably. The thickness of the separated region, 
however, begins to increase, as illustrated in figures 15 and 16. The remainder of the 
stream -function contours in figures 15(a) and (b) illustrate the growing thickness of the 
separated region, the increasing back flow, and the separation of flow eddies. At 
t - 2.23 a bubble begins to form on the trailing edge. This bubble continues to grow and 
is followed by the formation of another bubble at t = 2.53. The extent of both bubbles 
continues to increase until a single bubble is formed at t - 3. 13. The final • \}/ contours , 
given indicate that the single bubble has become much larger. The thickness of the s^-. 
aration region is roi^hly 1.25 times the airfoil thickness at this time, with the. forward'^ 
separation point at approximately 31 percent of the chord' 

The separation process may also be examined by viewing the series of upper - 
surface velocity profiles given in figure 16. The profiles given in figure 16(a) illustrate 
attached vqpper -surface flows shortly after the impulsive start. There is a noticeable 
increase in the boundary-layer thickness at t = 0.336. The remaihii^ parts of figure '16 
track the growii^ thickness of the separated region and the vq) stream movement of the 
separation point. 

Pressure and force coefficients at an early time are illustrated in figure 17(a). 

The friction drag constitutes more than 65 percent of the total drag at this time. As the 
boundary la]^r becomes thicker, the friction drag decreases rapidly and causes a cor- 
responding decrease in the total drag. The effects of the laminar floW separation are 
shown in figvire 17(b). The peculiar variation in the Cp distribution at the trailing 
edge is again due to inaccuracies in the body vorticity distribution. 
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Cambered Rock 

To show that the natural coordinate method could be used with arbitrary shaped 
bodies, the viscous flow about the cambered rock at a Reynolds number of 500 was 
develc^edi The contracted natural coordinate system used in the solution is given in 
fig^e 18. The same amplitude and decay factors given for the G&ttingen 625 airfoil 
Were used to create this system. Other summary data are presented in table 1. 

With a body such as the cambered rock, reliance must be placed on one's physical 

, I*- 

intuition in evaluating the resiiltii^ flow. For this reason an extended discussion of the 
flow past the cambered rock will not be given. Instead a significant number of xl/ and 
o» contours will be given. These are shown in figure 19. However, it is felt that some* 
remarks are appropriate. A glance at figure 18 indicates that the rock possesses sev- 
eral concave areas. Intuition would imply that flow stagnation areas should develop quite 
rapidly in these regions. This, in fact, does occur, as figure 19 indicates. In addition 
one would expect laminar flow separation and the consequent shedding of vorticity from'^^r 
the body. These events are also borne put by the contours. The ever-increasing size of 
the region of significant vorticity is quite apparent from the figures. Finally, velocity 
profiles and surface pressure distributions aire shown in figures 20 and 21. ■ .. . 

• i . Computer Time Requirements 

Numerical solutions to parabolic partial differential equations require extensive 
amounts of (%ital computer time. The total CPU times (UNIVAC 1106) used to generate 
the three solutions discussed in this study are documented in table 1. The average time 
required to converge each time step is also shown. (No attempt was made to quantify the 
effects of time -step size on the average times given.) . ' ' 

CONCLUSIONS 

The objective of this study was to develop methods to obtain numerical solutions of 
the two-dimensional, incompressible, time -dependent Navier -Stokes equations about 
arbitrary.bodies. The solutions followed the development of a general numerical curvi- 
linear. coordinate. transformation which produces a natural coordinate system having: a 
constant.- coordinate line coincident with each. boundary contour in the physical plane. 

Once the natural coordinates are developed for a given physical domain, the set of partial 
differential equations of interest may be transformed to the natural system and solved 
numerically in the transformed plane without regard to the geometry of the physical 
region. In effect the natural coordinate method eliminates all geometrical considera- 
tions from a given solution, as all physical regions have the same appearance in the 
transformed plane. The computer software utilized to generate the natural coordinates 
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is independent of the set of partial differential equations whose solution is to be carried 
out on the transformed plane. The partial differential equations governing potential and 
viscous flow differ drastically. However, for a given body geometry, the same natural 
system was used herein for both solutions. The second major advantage of using natu- 
ral coordinates is that the computer software generated to approximate the solution of a 
given set of partial differential equations is completely independent of the physical geom- 
etry of the problem. The same computer program was utilized to develop all the solu- 
tions for the wide variety of bodies discussed herein. Only the input varies with the 
body. Such a procedure obviously has significant ramifications in numerical mathemat- 
ics and all areas of physical science. 

Significant viscous -flow results were obtained for three different bodies, These 
included two general airfoil sections and one completely arbitrary body. The airfoil 
solutions developed computational wiggles near the sharp trailing edge at the start of the 
impulsive flow. The wiggles were generated as a result of large vorticity gradients 
which appeared in this region at the start. This disturbance was convected away from 
the body essentially undamped, but produced no significant disturbance neair the body at 
later times. The solutions show the formation and development of the boundary layer, 
laminar separation bubbles, and completely separated flow. Present results extend to a 
Reynolds number of 2000. Although the magnitude of the calculated force coefficients 
cannot be compared with experimental data, as none exist at this low Reynolds number, 
the time variation of these parameters agreed quite well with the flow pattern develop -r 
ment. The cambered-rock solution proved that the natural coordinate methods could be 
applied to very general bodies. The manner of this flow development also agreed with 
intuitive physical reasoning. 

There appears to be no basic barrier to higher Reynolds number solutions, since 
the means are at hand to contract the coordinate system about the body as much as 
desired. Preliminary runs at a Reynolds number of 10 000 have already been made. 

Work is also in progress on the solution for multiple airfoils with viscous flow. J ' 
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TABLE 1.- SUMMARY OF DATA FOR VISCOUS-FLOW SOLUTIONS 


Body 

(*) 

0, 

R 

Field 

Convergence criteria 

Time steps 

Average number 
of iterations 

Total 

Total 

deg 

size 


U) 

dxjy 

dn 

Initial 

Final 

Initial 

Final 

solution 

time 

CPU 
time, hr 

1 

15 

200 

3828 
(66 X 58) 

10-5 

10-5 

10-2 

0.01 

0.03 

280 

280. 

1.15 

7.i04 

2 

5 

2000 

4350 
(75 X 58) 

10-5 

10-5 

10-1 

.001 

.01 

.130 

220 

3.35 

32.40 

3 

5 

500 

3132 
(54 X 58) 

10-5 

10-5 

0.5 

.005 

.005 

55 

95 

1.30 

6.40 


Body code: 

1 Flapped K^m^-Trefftz airfoil 

2 Gottingen 625 airfoil 

3 Cambered rock • 


to 

<n 




























Transformed Plane 

Figure 2.- Field transformation - multiple bodies 





(b) Two-body region. 

Figure 3.- Computational grids - single and two-body regions 













LIFT COEFFICIENT 


2.538680 


□ NUMERICAL SOLUTION 


ANALYTIC S0LUTI0N(REF.6) 

DRAG COEFFICIENT = 0.004002 


LIF^ CCJEFFICIEN'!' ERROR = -0.00 IS 17 



Figure 6.” Analytic and numerical potential-flow results for flapped 

Karman-Tref ftz airfoil. 
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□ NUMERICAL SOLUTION 
ANALYTIC SOLUTION(REF. 6) 



(b) Streamlines . 
Figure 6.- Concluded 
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MOMENT COEFFICIENT = 


0.244996 
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(b) .Pressure distribution - fore body. 
Figure 10.- Continued. 
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(c) Pressure distribution - aft body. 
Figure . 10 Concluded . 
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Figure 12.- Velocity profiles for flapped Karman-Trefftz airfoil 


















t-0.658 


(a) t = 0.118, 0.336, and 0.658. 

Figure 16.- Upper-surface velocity profiles for Gottingen 625 airfoil 





t«1.83 


(b) t » 1.012, 1.53, and 1.83 
Figure 16.- Concluded. 
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PRESSURE DISTRIBUTION 
gStTINGEN 625 RIRFOIL 

RNGLE OF flTTflCK=5.00 • REYNOLDS . NUMBER=2000 

TIHE=3.33D0 
Cl»0. 415316 

Cn»0. 229930 Cn =0.153614 Cn =0.076316 













PRESSORE-'OISTRIBUTION 

CfiKBEREO'RQCK 

ANGLE OF flTT9CK=5.00 REYNOLDS NUMBER=500 

TIME=0.1000 
C, =-0.281309 

U. 


Cn=5. 439651 Cn =4.91361? Cp =0.526034 
^ D Dp Dp 

o 



(a) t = 0.1. 

Figure 21.- Pressure distribution for cambered rock. 




PRESSURE DISTRIBUTION 
CAMBERED ROCK 

ANGLE OF ATTACK=5.00 REYNOLDS NUMBER=500 

TIME=1.2000 
C|_=-0. 11B6146 

Ck= 0. 624591 Cn =0.512546 Cn=0. 112046 




